#include "matrices/matrix.h"
#include "comm/log.h"

#include <gsl/gsl_blas.h>
#include <gsl/gsl_linalg.h>
#include <json/json.h>


using namespace nymphaea;
using namespace nymphaea::matrices;

Matrix::Matrix(size_t row, size_t col):
	mpMatrix(NULL),
	mpInverse(NULL),
	mpLUDecomp(NULL),
	mpPerm(NULL),
	mDeterminant(0.0F),
	mRowSize(row),
	mColSize(col)
{
	mpMatrix = gsl_matrix_calloc(mRowSize, mColSize);
}

Matrix::Matrix(const Matrix& ins):
	mpMatrix(NULL),
	mpInverse(NULL),
	mpLUDecomp(NULL),
	mpPerm(NULL),
	mDeterminant(0.0F),
	mRowSize(ins.get_row_size()),
	mColSize(ins.get_col_size())
{
	mpMatrix = gsl_matrix_calloc(mRowSize, mColSize);
	gsl_matrix_memcpy(mpMatrix, ins.mpMatrix);
	
}

Matrix::Matrix(const gsl_matrix* const data):
	mpMatrix(NULL),
	mpInverse(NULL),
	mpLUDecomp(NULL),
	mpPerm(NULL),
	mDeterminant(0.0F),
	mRowSize(data->size1),
	mColSize(data->size2)
{
	mpMatrix = gsl_matrix_calloc(mRowSize, mColSize);
	gsl_matrix_memcpy(mpMatrix, data);
}

Matrix::~Matrix()
{
	reset_inverse_LUDecomp_deter_perm();
	gsl_matrix_free(mpMatrix);
	mpMatrix = NULL;
}

void Matrix::reset_inverse_LUDecomp_deter_perm()
{
	if(NULL != mpInverse){
		gsl_matrix_free(mpInverse);
		mpInverse = NULL;
	}

	if(NULL != mpLUDecomp){
		gsl_matrix_free(mpLUDecomp);
		mpLUDecomp = NULL;
	}
	
	if(NULL != mpPerm){
		gsl_permutation_free(mpPerm);
		mpLUDecomp = NULL;
	}
	mDeterminant = 0.0F;
}

size_t Matrix::get_row_size() const
{
	return mRowSize;
}

size_t Matrix::get_col_size() const
{
	return mColSize;
}

Matrix Matrix::get_row_by_idx(size_t idx) const
{
	if(idx > get_row_size()){
		return Matrix(1, get_col_size());
	}
	gsl_vector* tmp_vector = gsl_vector_calloc(get_col_size());
	gsl_matrix* tmp_matrix = gsl_matrix_calloc(1, get_col_size());

	gsl_matrix_get_row(tmp_vector, mpMatrix, idx);
	gsl_matrix_set_row(tmp_matrix, 0, tmp_vector);

	Matrix ret(tmp_matrix);

	gsl_vector_free(tmp_vector);
	gsl_matrix_free(tmp_matrix);

	return ret;
}


Matrix Matrix::get_col_by_idx(size_t idx) const
{
	if(idx > get_col_size()){
		return Matrix(get_row_size(), 1);
	}
	gsl_vector* tmp_vector = gsl_vector_calloc(get_row_size());
	gsl_matrix* tmp_matrix = gsl_matrix_calloc(get_row_size(), 1);

	gsl_matrix_get_col(tmp_vector, mpMatrix, idx);
	gsl_matrix_set_col(tmp_matrix, 0, tmp_vector);

	Matrix ret(tmp_matrix);

	gsl_vector_free(tmp_vector);
	gsl_matrix_free(tmp_matrix);

	return ret;
}

double Matrix::get_element_by_idx(size_t r_idx, size_t c_idx) const
{
	return gsl_matrix_get(mpMatrix, r_idx, c_idx);
}

void Matrix::set_element_by_idx(size_t r_idx, size_t c_idx, double val)
{
	reset_inverse_LUDecomp_deter_perm();
	return gsl_matrix_set(mpMatrix, r_idx, c_idx, val);
}

void Matrix::set_all(double val)
{
	reset_inverse_LUDecomp_deter_perm();
	gsl_matrix_set_all(mpMatrix, val);
}


Matrix Matrix::get_transposed() const
{
	gsl_matrix* transposed = gsl_matrix_calloc(get_col_size(), get_row_size());
	for(size_t i = 0; i < get_row_size(); ++i){
		for(size_t j = 0; j < get_col_size(); ++j){
			gsl_matrix_set(transposed, j, i, gsl_matrix_get(mpMatrix, i, j));
		}
	}
	return Matrix(transposed);
}

Matrix Matrix::get_inverse() 
{
	if(get_row_size() != get_col_size()){
		// Not square matrix
		return Matrix(1,1);
	}
	if(NULL == mpInverse){
		// Cal inverse matrix once

		cal_LUDecomp_deter_perm();
		mpInverse = gsl_matrix_calloc(get_row_size(), get_col_size());
		gsl_linalg_LU_invert(mpLUDecomp, mpPerm, mpInverse);
	}
	return Matrix(mpInverse);
}

Matrix Matrix::get_lu_decomp()
{
	cal_LUDecomp_deter_perm();
	return Matrix(mpLUDecomp);
}

double Matrix::get_determinant()
{
	cal_LUDecomp_deter_perm();
	return mDeterminant;
}

double Matrix::sum() const
{
	double sum = 0.0F;
	for(size_t i = 0; i<get_row_size(); ++i){
		for(size_t j = 0; j<get_col_size(); ++j){
			sum += get_element_by_idx(i, j);
		}
	}
	return sum;
}


Matrix Matrix::operator*(const Matrix& other) const
{
	Matrix res(get_row_size(), other.get_col_size());
	gsl_blas_dgemm(CblasNoTrans, CblasNoTrans, 1.0F, (*this).mpMatrix, other.mpMatrix, 0.0F, res.mpMatrix);
	return res;
}

Matrix Matrix::operator/(const Matrix& other) const
{
	return Matrix(get_row_size(), other.get_col_size());
}

Matrix Matrix::operator+(const Matrix& other) const
{
	Matrix res(*this);
	gsl_matrix_add(res.mpMatrix, other.mpMatrix);
	return res;
}

Matrix Matrix::operator-(const Matrix& other) const
{
	Matrix res(*this);
	gsl_matrix_sub(res.mpMatrix, other.mpMatrix);
	return res;
}

bool Matrix::operator==(const Matrix& other) const
{
	if(get_row_size() != other.get_row_size()
			|| get_col_size() != other.get_col_size()){
		return false;
	}

	for(size_t i = 0; i < get_row_size(); ++i){
		for(size_t j = 0; j < get_col_size(); ++j){
			if(get_element_by_idx(i, j) != other.get_element_by_idx(i, j)){
				return false;
			}
		}
	}
	return true;
}
std::vector<double> Matrix::row_to_vector(size_t idx) const
{
	return std::vector<double>();
}

std::vector<double> Matrix::col_to_vector(size_t idx) const
{
	return std::vector<double>();
}

void Matrix::cal_LUDecomp_deter_perm()
{
	if(get_row_size() != get_col_size()){
		// Not square matrix
		return;
	}
	if(NULL != mpLUDecomp){
		// Cal Left Upper Decomp once
		return;
	}

	int s = 0;
	mpPerm = gsl_permutation_calloc(get_row_size());
	mpLUDecomp = gsl_matrix_calloc(get_row_size(), get_col_size());
	gsl_matrix_memcpy(mpLUDecomp, mpMatrix);
	gsl_linalg_LU_decomp(mpLUDecomp, mpPerm, &s);
	mDeterminant = gsl_linalg_LU_det(mpLUDecomp, s);
}

const std::string Matrix::dump_to_bytes() const
{
	Json::Value root;
	Json::StyledWriter writer;
	root["type"] = "dense";
	root["row"] = static_cast<unsigned int>(get_row_size());
	root["col"] = static_cast<unsigned int>(get_col_size());

	for(size_t i = 0; i<get_row_size(); ++i){
		for(size_t j = 0; j<get_col_size(); ++j){
			root["data"][i][j] = get_element_by_idx(i, j);
		}
	}

	return writer.write(root);
}

int Matrix::load_from_bytes(const std::string& data)
{
	Json::Reader reader;
	Json::Value root;

	if(!reader.parse(data, root)){
		LOG("Json parsing failed");
		// error log
		return -1;
	}
	size_t row = root.get("row", 0).asUInt();
	size_t col = root.get("col", 0).asUInt();

	mRowSize = row;
	mColSize = col;
	gsl_matrix_free(mpMatrix);
	mpMatrix = gsl_matrix_calloc(row, col);
	reset_inverse_LUDecomp_deter_perm();

	if(0 != row && 0 != col && root.get("type", "undefined") == "dense"){
		Json::Value data(root["data"]);
		if(row != data.size()){
			LOG("Row size not matched");
			return -1;
		}else{
			for(size_t i = 0; i<row; ++i){
				if(col != data[i].size()){
					// reset matrix to zero
					set_all(0.0F);
					LOG("Col size not matched");
				}else{
					for(size_t j = 0; j<col; ++j){
						set_element_by_idx(i, j, data[i][j].asDouble());
					}
				}
			}
		}
	}else{
		LOG("Matrix type is undefined.");
		return -1;
	}

	return 0;
}

